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We construct an efficient zero-temperature semi-local density functional to dynamically simulate 
an electron bubble passing through superfluid ''He under various pressures and electric fields up to 
nanosecond timescale. Our simulated drift velocity can be quantitatively compared to experiments 
particularly when pressure approaches zero. We find that the high-speed bubble experiences re- 
markable expansion and deformation before vortex nucleation occurs. Accompanied by vortex-ring 
shedding, drastic surface vibration is generated leading to intense phonon radiation into the liquid. 
The amount of energy dissipated by these phonons is found to be greater than the amount carried 
away solely by the vortex rings. These results may enrich our understanding about the vortex 
nucleation induced energy dissipation in this fascinating system. 
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I. INTRODUCTION 



An electron injected into liquid helium forms a bubble 
due to the Pauli exclusion between an excess electron 
and helium atoms. ^'^ The dissipation mechanisms of a 
moving electron bubble in superfluid '*He has attracted 
considerable research interests for many yearsi^""— Above 
1 K, the moving bubble experiences a drag force from 
collisions with thermally excited phonons and rotonsi^ 
Below 1 K, this drag force becomes very small, and even 
a weak electric field can accelerate the bubble to a high 
speed until some new dissipation mechanisms set in. In 
the high pressure regime, roton emission dominates the 
dissipation owing to the relatively low Landau velocity 
Vi^ and the high vortex-ring nucleation critical velocity 
Vc^^ In the low pressure regime, vortex-ring nucleation 
plays the key role. A vortex ring can attach to the bubble 
to form a bubble-ring complex if the electric field is not 
overly strong; otherwise, successive vortex rings can be 
shed away from the bubble surface.—""— 

Pioneering simulations using the Gross-Pitaevskii 
equation (GPE) have demonstrated the above vortex- 
ring nucleation, trapping, and shedding scenario J^ii^ 
But it is known to be difficult for the efficient local 
GPE to reproduce helium properties, ^^ In contrast, ac- 
curate nonlocal density functional (NLDF) theories have 
been extensively applied to quasiparticles, vortices, and 
ions-related problems in liquid helium J^ii^ But they usu- 
ally require prohibitively high computational cost for dy- 
namic simulations. In this paper, we first introduce our 
well constructed semi-local density functional (SLDF), 
associated with an optimized numerical scheme, that can 
reconcile both the physical accuracy and the computa- 
tional efficiency in its applicable regimej ^^d^ We then 
present our dynamic simulation for a moving electron 
bubble in pure superfluid ^He under low pressures and 
strong electric fields, where successive vortex-ring shed- 
ding is indeed observed. However, we shall point out 



that although vortex nucleation does trigger the dissipa- 
tion, the major part of energy loss may come from its 
induced phonon radiation via surface vibration, rather 
than purely the shed-away vortex rings. 



II. THEORETICAL FORMULATION 

We formulate our problem in the framework of a zero- 
temperature SLDF theory. The system free energy den- 
sity Q consists of the helium part, the electron part, and 
the helium-electron interaction part. 
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The helium part Ghc takes the form of 
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where ip is the macroscopic helium wavefunction and g = 
I^P is the local helium number density. The electron part 
Gc takes the form of 



2m, 



\V(t>f - eSzr], 
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where is the single electron wavefunction and 77 = 10 ^ is 
the local electron number density. The jVV'l^ and |V0p 
terms above are the helium and electron kinetic energy 
densities, with ttiho and m^ being the helium and electron 
masses, respectively, fj, is the helium chemical potential 
controlled by a given pressure p, which fixes the helium 
number density p in bulk. £ is the applied electric field, 
which drives the electron to move along z direction. 52, 
.93, 34 and /12, ^3, ^4 sue all fitting parameters to be 



explained below. The helium-electron interaction part 
Ghc-c takes the form of a contact collision 



^Hc-c = flQTli 



(4) 



in which /i is a fitting parameter chosen so as to produce 
a 1 A scattering length and hence a 1.0456 eV potential 
barrier for an electron to tunnel into homogeneous helium 
at zero temperature and zero pressure.^^ 

The polynomial function of g in the second line of 
Eq. ^ is the helium-helium interaction energy density in 
the form of two-particle, three-particle and four-particle 
contact collisions. In the homogeneous case, it gives the 
helium internal energy density 



from which one can derive the equation of state 



P[p\ = 2^2^^ 
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the chemical potential 

^^[p\ = 92P + 9zP^ +94P^ 
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and the sound velocity 

c[p] = V(.92/3 + 2g3p2 + 3g4P^) /mHc- (8) 

By choosing (72, 53 and 54 according to the well-known 
Orsay-Trento density functional, "'^^ the experimentally 
measured above quantities at zero temperature can be 
very well produced. 

The mixture polynomial function of g and |Vgp in 
the third line of Eq. ([2]) contains our introduced semi- 
local interactions also up to four-particle collisions. They 
are the lowest-order gradient-expansion corrections to the 
local interactions and are essential for incorporating the 
helium surface tension into our theory 
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By choosing /12, ^3 and ^,4 appropriately, the experimen- 
tally measured surface tension at zero temperature and 
zero pressure can be produced. In addition, Eq. ^ is 
a generalized formula to allow not only zero pressure, 
when the liquid is in equilibrium with vacuum, but also 
a positive pressure, when the liquid is in contact with an 
impenetrable wall. In both situations, the liquid density 
drops from its bulk value p to within a thin interfacial 



layer. a[p\ measures the energy change per unit surface 
area due to this density bending, in comparison with the 
internal energy per unit area held by the same number 
of particles in bulkfil. Such a generalized definition on 
the surface tension turns out to be useful for the elec- 
tron bubble problem, where the bubble boundary pushes 
away the liquid basically like an impenetrable wall, and 
so gives rise to a nontrivial surface energy under any fi- 
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TABLE I: The chosen fitting parameters 


of our SLDF. 
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FIG. 1: (color online). The helium properties in our SLDF. 
(a) The equation of state, (b) The surface tension versus pres- 
sure, (c) The excitation spectrum under the spinodal pressure 
—9.5 bar, the saturation pressure bar, and the melting pres- 
sure 25 bar. (d) The vortex-line energy per unit length versus 
pressure with a 20 A cutofT radial distance. 

Table U lists all the chosen fitting parameters of our 
SLDF. Figure [Ija) shows the equation of state. At 
p = bar, po = 0.021836 A'^, /zq = -7.1500 K, and 
Co = 237.70 m s^^. Figure[IJb) shows the surface tension 



versus pressure. At p = bar, (Tq = 0.37554 erg cm~^ 
with a calculated 5.8 A surface thickness (defined by 
10% to 90% helium density) consistent with experiments 
and other theories J^iiiiiS 

The excitation spectrum in our theory can be found to 
be of the Bogoliubov type 



Putting them into the hollow-core modeliS 
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whose pressure dependence is enclosed in the sound veloc- 
ity c[p] from the local interactions, and the dimensionless 
function b[p] from the semi-local interactions. Compared 
with the efficient local GPE, although our SLDF can 
indeed produce more realistic helium bulk and surface 
properties, and meanwhile maintain the same computa- 
tional efficiency, it still cannot incorporate the backfiow 
effect and the roton excitation unless some nonlocal in- 
teractions are introducedJ^ii^ Figurelljc) shows the exci- 
tation spectrum under three typical pressures: the spin- 
odal pressure —9.5 bar, the saturation pressure bar, 
and the melting pressure 25 bar. The primary reason 
for employing three semi- local parameters, /12, /13 and 
/14, in our SLDF instead of a single one, as that done in 
the well-known Stringari-Treiner density functional, "'^^ is 
to remove some pathological behavior of the excitation 
spectrum. A single semi-local parameter being used to 
fit the surface tension at zero pressure will result in too 
large of a b[p] and hence unrealistically high excitation 
energies even at just moderate momentumsJ^ This will 
make the liquid overrigid against density fiuctuation and 
may bring on artifacts in dynamic simulations. With the 
three semi-local parameters, we can keep b[p] sufficiently 
small, and also let it slowly decrease with pressure to 
imitate the tendency of roton gap with pressure in real 
helium. 

We can also find in our theory the rectilinear vortex- 
line energy per unit length, by calculating the radial num- 
ber density distribution g{r) under a given pressure. 
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where g drops from its bulk value p to as r — >■ 0, and 
g' = dg/Ar. The first term in the first line contributes 
to the kinetic energy Kn^c/L due to one unit of quan- 
tum circulation. All the other terms contribute to the 
potential energy Unne/L due to the density bending. ^° 
The critical negative pressure for the free expansion of a 
vortex core is found to be close to —7 bar. Figure [TJd) 
shows the vortex-line energy per unit length versus pres- 
sure with a d = 20 A cutoff radial distance. At p = bar, 
we find ifii„c/£ = 2.94 K A^^ and Unnc/L = 1.57 K A^^. 
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we can get the core size a — 0.581 A and the core param- 
eter S = 1.893, which imply a rather rigid core. 



III. DYNAMIC SIMULATION 

We perform the dynamic simulation in the experimen- 
tal low pressure regime p = ~ 10 bar and strong electric 
field regime £ = 1 ~ 30 MV m~^. Under these con- 
ditions, a high-speed electron bubble can nucleate vor- 
tex rings. As is well-known, there are two classic mod- 
els on the formation of a vortex ring and bubble-ring 
complex.—^— One is the encircling-ring (quantum tran- 
sition) model, in which a full ring axisymmetrically ap- 
pears on the bubble equator first, then moves sideways 
and captures the bubble on its core. The other one is 
the pinned-ring (peeling) model, in which a small proto- 
ring non-axisymmetrically appears on the bubble equa- 
tor first, then grows up and retains the bubble on its 
core. Experimental evidence suggests that the encircling- 
ring model gives a better description at lower tempera- 
tures in isotopically purified ^He, whereas the pinned- 
ring model gives a better description at higher temper- 
atures in natural He, with the existence of thermal ro- 
tons or '^He impurities. ^^^ For an already formed bubble- 
ring complex, the strength of the electric field deter- 
mines its subsequent behavior. In a field of the order 
of 10"^ ^ lO** V m~^, the ring can keep trapping the 
bubble and may grow to micron size, which significantly 
hinders the bubble motioni^ But in a field of the or- 
der of 10^ ~ 10^ V m~^, owing to the vanishing escape 
barrier, the ring should separate from the bubble dur- 
ing several picoseconds with its size almost unchanged. 
The bubble can then move more freely, as well as suc- 
cessively shed away vortex rings.^ Since our simulation 
is for pure '*He at zero temperature, the axisymmet- 
ric vortex-ring nucleation in the encircling-ring model 
is more suitable. Furthermore, since our electric fields 
lie in the strong regime, the transient non-axisymmetric 
capture-and-escape processes should have little influence 
on the longtime physics. Thus we reduce the original 
three-dimensional problem into two-dimensional by tak- 
ing cylindrical symmetry with transverse coordinate r 
and longitudinal coordinate z. This simplification makes 
the numerical integration up to nanosecond timescale 
computationally practical. 

The coupled equations of motion for the macroscopic 
helium wavefunction ip and the single electron wavefunc- 



tion 6 can be derived from our SLDF as 
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The fourth-order finite-difference method in space and 
the fourth-order Runge-Kutta method in time are used 
to do the numerical integration. The computation grid 
on the r-z plane is 500 x 1000 with the space step of 
0.2 A. So the physical space is a pipe of 200 A in diameter 
and 200 A in length. We let it move along to keep the 
electron in the center. Systematic numerical tests have 
been done to ensure that 0.2 A spatial resolution provides 
enough precision for our problem, even when the vortex 
core structure is involved. 

Due to the big difference between the helium and elec- 
tron intrinsic timescales of the order of mno/'Tie ~ 7300, 
we adopt the adiabatic approximation. For every instan- 
taneous helium configuration, we first find the electron 
ground state by evolving the electron with 50 imaginary 
time steps of 0.001 fs, then develop helium with 1 real 
time step of 0.001 ps. To avoid the outgoing sound waves 
from the central region being reflected on the grid bound- 
ary, we set a space-dependent damping coefficient by re- 
placing i with ri{s) -|- i in the helium equation of motion, 
where 
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{s = v/r-2+z2). (15) 



We simply choose s = 80 A, w = 1.0 A, and fj = 0.5. 
The sound waves can travel freely through the undamped 
region {s < s) but are quickly attenuated in the damped 
region {s > s). Physically, this looks as if there is a strong 
zero-temperature heat sink surrounding the system. 

Figure [2] shows the calculated radial number density 
profiles of helium g{r)/po and electron r]{r)/r]o for a 
spherical static electron bubble under various pressures. 
Here rjo = t^/'^R^ with _Ro = 20 A being merely a refer- 
ence radius. The static bubble radius Rh (defined at 50% 
helium density) can be found to be 18.52 A at p = bar 
and squeezed down to 13.50 A at p = 10 bar, consis- 
tent with experiments and other theories.-!'- We set these 
ground states as the initial states to start the time evo- 
lution. For every single run with a specified pressure and 
electric field, 8-core parallel computation at 2.83 GHz was 
performed. Multiple runs with different pressures and 
electric fields were conducted simultaneously on our com- 
puting cluster which possesses hundreds of cores. Typ- 
ically, it took about one month to reach 10 ns physical 
time. 
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FIG. 2: (color online). The radial number density profiles 
of helium and electron for a spherical static electron bubble 
under various pressures. 



Figure [3] shows the snapshots of the helium number 
density profile Q{r, z; t)/po around the moving electron 
bubble at a sequence of times under p = bar and 
£ = IQ MV m~^. The instantaneous velocity v at each 
time is approximated by linearly fitting the displacement 
of the electron wavefunction-maximum over a 20 ps time 
interval. At i = ps, the bubble is on its spherical ground 
state of radius i?b = 18.52 A, when the driving field is 
just switched on. The initial acceleration of the bubble 
has the value that is expected based on the applied force 
e£ and the hydrodynamic mass 27ri?j^r7iHcPo/3. Up until 
t = 40 ps, V is already as high as 30.6 m s^^, whereas the 
bubble shape is only slightly squeezed and the liquid in 
front is slightly compressed. As the bubble velocity in- 
creases further, its volume keeps expanding and its shape 
keeps deforming into an oblate spheroid due to the nega- 
tive pressure on the bubble waist caused by the Bernoulli 
effect.— The acceleration decreases remarkably because 
of the increase in the effective mass so that v approaches 
a constant value about 50.6 m s~^ and the flow pattern 
exhibits quite well fore-and-aft symmetry as shown at 
t — 450 ps. After this stage, a narrow edge girdling the 
bubble waist gradually sticks out signifying the onset of 
vortex-ring nucleation. This adds a considerable inertia 
to the bubble motion, making v decrease to 44.9 m s~^ 
at t = 850 ps, when the bubble attains its largest equa- 
torial radius about 44 A and its sharpest edge radius of 
curvature about 4 A. We may call the bubble velocity for 
the first appearance of this characteristic behavior as the 
vortex-ring nucleation critical velocity v^ — 44.9 m s~^ 
under the prescribed conditions. There is no dissipa- 
tion before, since the electric work mainly transforms to 
the fluid kinetic energy. But around this time the liquid 
density starts to fluctuate accompanied by some phonon 
generation, presumably because the local flow near the 
bubble waist becomes supersonic and keeps hitting the 
nucleated vortex core i"^^i"^^ At t = 950 ps, a vortex ring 
of radius i?ring ~ 36 A is detaching from the bubble sur- 
face, while V is 43.8 m s~^ and continues to slow down 
due to the attraction from the vortex ring. The phonon 
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FIG. 3: (color online). The snapshots of the helium number density profile around the moving electron bubble at a sequence 
of times under zero pressure and 10 MV m~^ electric field.^ The bubble travels from left to right along 2-direction (the 
longitudinal direction). Each image corresponds to a physical size of 100 A x 100 A and is reflected up and down in r-direction 
(the transverse direction) justified by the imposed cylindrical symmetry. The instantaneous velocity is given for each time. 
The open and solid arrows indicate the movements of vortex ring and solitary ripplon relative to the bubble as mentioned in 
the text. 



generation appears more intense, particularly peaked at 
a wavelength of the order of 10 A. The unexpected large 
ring radius, almost twice as big as the static bubble ra- 
dius, is a unique observation of our dynamic SLDF sim- 
ulation. 

In the later stage, the vortex ring falls behind, leav- 
ing a drastic surface vibration and velocity oscillation 
to the bubble in the next hundreds of picoseconds. At 
t — 970 ps, an annular surface protrusion like a soli- 
tary ripplon is formed and sliding backwards, while v 
is dragged down to 31.5 m s~^. This solitary ripplon 
then shrinks on the bubble tail at i = 990 ps, and radi- 
ates a large number of phonons backwards into the liq- 
uid as shown a.t t — 1030 ps. This provides a strong 
forward impulse to speed up the bubble so that v can 
even overshoot to 60.5 m s~^ at i = 1070 ps. During 
this time, the solitary ripplon surviving from the last 
collision bounces back and propagates forwards. Around 
t = 1080 ~ 1090 ps, it shrinks on the bubble head and de- 
posits a large number of phonons forwards into the liquid. 
This provides a strong backward impulse decelerating the 
bubble velocity to as low as 27.7 m s~^ at i = 1170 ps. 
Afterwards, some remnant surface vibration and phonon 
radiation keep happening back and forth, while the bub- 
ble velocity oscillatorily increases until a second vortex 
ring completes its nucleation. Then all the processes re- 
peat. 



IV. DISCUSSION 

Figure |4] shows the time evolution of the moving elec- 
tron bubble volume Vb(i)/(47ri?j^/3) under p = bar 
and f = 10 MV m"-'^, obtained from the same evolu- 
tion shown in Fig. [3] but extended to about 5 ns physical 
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FIG. 4: (color online). The time evolution of the moving elec- 
tron bubble volume under zero pressure and 10 MV m~^ elec- 
tric field. The dotted lines mark the moments of the vortex- 
ring shedding events. 



time. Here 47rJ?^/3 with i?b = 18.52 A is the spherical 
static bubble volume at i = ps. As can be seen, be- 
fore the dissipation processes arise at t « 850 ps, Vh{t) 
keeps growing to about 3.6 times the static volume. After 
that, Vh{t) oscillates violently between 2.0 and 4.0 (aver- 
age about 3.0) times the static volume. It approximately 
implicates the hydrodynamic mass variation, and coin- 
cides with the bubble velocity variation described above. 
The dotted lines mark the moments of the vortex-ring 
shedding events. During t = 800 ~ 5000 ps, there are six 
such events in total. Although they do not exhibit a regu- 
lar periodicity, we may still roughly estimate the average 
time interval between two successive events r « 700 ps. 
Each time when a ring leaves the bubble, there follows 
a huge volume change immediately and perhaps a series 
of small volume changes later on. This can be related to 



the complex surface motion displayed in Fig. |3l 
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FIG. 5: (color online). The drift velocity versus pressure 
under 2.6, 10, 20, 30 MV ni~ electric fields from our sim- 
ulation in comparison with the low pressure part of experi- 
ment by Nancolas et. al. at 0.3 K under 2.6 MV m~^. The 
dashed curve gives the Landau velocity for roton emission. 
The dashed-dotted curve gives the theoretically extrapolated 
drift velocity from 25 bar under 2.6 MV m~^ if roton emission 
was the only dissipation mechanism.— 

Figure [5] shows the drift velocity v versus pressure un- 
der £ = 2.6, 10, 20, 30 MV m~^ from our simulation in 
comparison with the low pressure part of experiment by 
Nancolas et. al. at 0.3 K under 2.6 MV m'^^ Our v 
is calculated by linearly fitting the displacement of the 
electron wavefunction-maximum over 2 ns after the bub- 
ble motion has become strongly dissipative. Although 
the instantaneous velocity v mentioned above undergoes 
large oscillation, v is nearly constant depending only on 
the pressure and electric field for a long time. According 
to Nancolas et. al. , if roton emission was the only dissipa- 
tion mechanism, then in the limit £ — > V m^^, v should 
follow the dashed curve of Landau velocity Wl; whereas 
for £ = 2.6 MV m~^, v should follow the dashed-dotted 
curve based on the well-known relation v = Vi^+A£^'^ , in 
which A takes its value at 25 bar with the pressure depen- 
dence ignored. However, the experimentally measured 
V shows a significant drop from this curve, indicating 
an additional dissipation mechanism in this low pressure 
regime believed to be vortex-ring nucleation. Our work 
supports this general scenario especially for p < 3 bar 
where the simulation and experimental results quanti- 
tatively agree with each other. This suggests that our 
SLDF can be an efficient and reliable tool to study inter- 
esting dynamic problems, particularly under the natural 
experimental condition p — > bar close to zero tempera- 
ture. The lack of roton excitation in our SLDF may not 
be problematic for two reasons that essentially prohibit 
roton emission in this range. First, the Landau velocity 
^L ~ 58 m s~^ is substantially higher than the vortex- 



ring nucleation critical velocity Wc ~ 45 m s~^. Sec- 
ond, there is experimental evidence that A is divergent 
as p — > 3 bar from above.— For p > 3 bar, our simulation 
still exhibits a qualitatively correct behavior in terms of 
the increasing v with increasing p commonly understood 
as the inverse proportionality of Vc with respect to the 
bubble size.—"— The relatively large deviation between 
simulation and experiment is presumably caused by the 
lack of roton excitation in our SLDF, whereas in reality 
vortex-ring nucleation and roton emission coexist in this 
pressure range. 

Although in Fig. |3]we have presented that vortex nu- 
cleation is indeed the origin of dissipation, it is nontrivial 
to investigate how the energy is actually taken away from 
the system. While this is normally attributed to purely 
the frequent vortex-ring shedding, our analysis suggests 
otherwise. The energy associated with a large single vor- 
tex ring at zero pressure can be estimated bj»i^ 
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Based on our preceding calculation about the vortex 
properties in our density functional, when i?ring ~ 36 A 
as shown in Fig. [3l we have E'ring ~ 1147 K. As a re- 
sult, the average rate of energy loss to the vortex rings is 
^ring/rsa 1147K/700ps = 2.3xl0-"W. On the other 
hand, the average rate of energy input to the system from 
the electric work is e£v « 10 MeV m~ x 48 m s"'^ = 
7.7 X 10~^^ W. Strikingly, the released vortex rings only 
take away about 30% of the total input energy. Accord- 
ing to the processes illustrated in Fig. |3] and discussed 
above, the remaining energy loss can be primarily at- 
tributed to the phonon radiation through surface vibra- 
tion persistently occurring in the time interval between 
successive vortex-ring shedding events. Under a fixed 
electric field with an increasing pressure, we even notice 
a tendency for such a dissipation channel to become in- 
creasingly dominant. 



V. CONCLUSION 

We have carried out dynamic simulation for a moving 
electron bubble in pure superfiuid "^He under low pres- 
sures and strong electric fields, using a zero-temperature 
semi-local density functional theory that embodies re- 
alistic ""^He bulk and surface properties. It visualizes 
the underlying microscopic processes from picosecond to 
nanosecond timescale, which cannot yet be directly ob- 
served through any experimental apparatus today. Not 
only have we confirmed that theory and experiment show 
good agreement on the bubble mobility approaching zero 
pressure, we have also discovered that the pronounced 
surface vibration and phonon radiation as a result of vor- 
tex nucleation dissipate more energy than the shed-away 
vortex rings alone. These results may enrich our under- 
standing about the vortex nucleation induced energy dis- 



sipation in this fascinating system and so warrant further 
experimental and theoretical studies. 
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